%% test script
%
%%
clear all
clc
%%
img=imread('lena.jpg');
imgsz=size(img)

figure(1);
imshow(img);

row = imgsz(1);
col = imgsz(2);

ut=1

%Mn=[img(:,:,1)-repmat(mean(img(:,:,1),2),1,col);img(:,:,2)-repmat(mean(img(:,:,2),2),1,col);img(:,:,3)-repmat(mean(img(:,:,3),2),1,col)];
Mean1=repmat(mean(img(:,:,1),2),1,col);
Mean2=repmat(mean(img(:,:,2),2),1,col);
Mean3=repmat(mean(img(:,:,3),2),1,col);

Mn1=img(:,:,1)-Mean1;
Mn2=img(:,:,2)-Mean2;
Mn3=img(:,:,3)-Mean3;

[U1,S1,V1]=svd(Mn1);
[U2,S2,V2]=svd(Mn2);
[U3,S3,V3]=svd(Mn3);

B1=zeros(row,col);
B2=zeros(row,col);
B3=zeros(row,col);

sumSi1=0;
sumSi2=0;
sumSi3=0;

SumS1=sum(sum(S1));
SumS2=sum(sum(S2));
SumS3=sum(sum(S3));

for i=1:size(S1,1)

	sumSi1=S1(i,i)+sumSi1;	
	if sumSi1 >=ut*SumS1
                break;
        end
end
diagnum1=i
B1=U1(:,1:i)*S1(1:i,1:i)*V1(:,1:i)';

for i=1:size(S2,1)

        sumSi2=S2(i,i)+sumSi2;
        if sumSi2 >=ut*SumS2
                break;
        end
end
diagnum2=i
B2=U2(:,1:i)*S2(1:i,1:i)*V2(:,1:i)';

for i=1:size(S3,1)

        sumSi3=S3(i,i)+sumSi3;
        if sumSi3 >=ut*SumS3
                break;
        end
end
diagnum3=i
B3=U3(:,1:i)*S3(1:i,1:i)*V3(:,1:i)';


imgnew(:,:,1)=floor(B1+Mean1);
imgnew(:,:,2)=floor(B2+Mean2);
imgnew(:,:,3)=floor(B3+Mean3);

%imgnew(imgnew>255)=255;
imgnew=uint8(imgnew);

figure(2)
imshow(imgnew)
Sz=img-imgnew;
